function pop_sac()
%
% compute the SACs across trials for cells in a population.
% 16 Jan 06 P. Manis
%
% to do this, we rearrange the data to feed it into the SAC routine.
% 
global CONTROL

sf=getmainselection; % get the list.
pflag = getplotflag;

twin = 100;
binw = 0.02;
binw2 = 0.1;
start1 = 100;
dur1 = 800;

X={};
Y={};
ntrial =  3;
trial = 1;
for i = 1:length(sf)
    sfi = sf(i);
    if(isfield(CONTROL(sfi).spike, 'latency'))
        X = CONTROL(sfi).spike.latency; % data from all 100 trials
        Y{i} = X{trial};
    else
        Y{i}=[];
    end;
    
end;
[y, yh1, hx1, mr1] = sac(Y, twin, binw2, start1, dur1); % different bin width for single trace.
printrasterplotfile('rasterplot1.txt', 1, Y); % subroutine at the end of this source file

y2s=[];
yh2s = [];
hx2s = [];
mr2s = [];
for trial = 1:ntrial;
    Y={};
    for i = 1:length(sf)
        sfi = sf(i);
        if(isfield(CONTROL(sfi).spike, 'latency'))
            X = CONTROL(sfi).spike.latency; % data from all 100 trials
            Y{i} = X{trial};
        else
            Y{i}=[];
        end;
    end;
    [y2, yh2, hx2, mr2] = sac(Y, twin, binw, start1, dur1);
    %    y2s(trial,:) = y2;
    yh2s(trial,:) = yh2;
    hx2s(trial,:) = hx2;
    mr2s(trial) = mr2;
    if(trial == 2)
        printrasterplotfile('rasterplot2.txt', 2, Y); % subroutine at the end of this source file
    end;

end;
% y2 = mean(y2s);
yh2 = mean(yh2s);
hx2 = mean(hx2s);
mr2 = mean(mr2s);



% perform gaussian fits on the histograms to identify the peaks and measure
% them. Using Molitor's routines (mrqfit).
% [...] = MRQFIT(F, P, X, Y, SIG, VP, LB, UB, IMAX, TOL)
% F = 'gaussian', P = [A0 A1 M1 S1 A2 M2 S2 ... AN MN SN]
% SIG is sigma Y (def = 1); VP is vary array per parameter [0 or 1];
% LB, UB are upper and lower bounds.
%
hx1c = hx1+0.5*(hx1(2)-hx1(1));
hx2c = hx2+0.5*(hx2(2)-hx2(1));

gpar1 = [1, 5, 0, 1]; % single gaussian centered on 0
scf = 1;
gulim1 = [0, 100, 10*scf, 20*scf];
gllim1 = [0, 0, 0.0, 0.01];
gvar1 =  [0 1 0 1];
nitermax = 100;
gpar2 = [gpar1 5 10*scf 1*scf];
gulim2 = [gulim1 100 50*scf 20*scf];
gllim2 = [gllim1 0 2*scf 0.1*scf];
gvar2 = [gvar1 1 1 1];

[fp11, chisq11, niter11, fitc11, err11, dep11] = mrqfit('gaussian', gpar1, hx1c*scf, yh1, [], gvar1, gllim1, gulim1, nitermax, []);
[fp21, chisq21, niter21, fitc21, err21, dep21] = mrqfit('gaussian', gpar2, hx1c*scf, yh1, [], gvar2, gllim2, gulim2, nitermax, []);

[fp12, chisq12, niter12, fitc12, err12, dep12] = mrqfit('gaussian', gpar1, hx2c*scf, yh2, [], gvar1, gllim1, gulim1, nitermax, []);
[fp22, chisq22, niter22, fitc22, err22, dep22] = mrqfit('gaussian', gpar2, hx2c*scf, yh2, [], gvar2, gllim2, gulim2, nitermax, []);


yg1 = gaussfunc(hx1c, fp11);
yg2 = gaussfunc(hx2c, fp12);
fwhmfac = 2*sqrt(2*log(2)); % note - log is ln.
fwhm1 = fwhmfac*fp11(4); % full width at half maximal height.
fwhm2 = fwhmfac*fp12(4);


hf = findobj('tag', 'SAC2');
if(isempty(hf))
    hf = figure;
    set(hf, 'Tag', 'SAC2');
end;
% generate the results structure: SAC
SAC.hx1 = hx1; % save the histograms.
SAC.hy1 = yh1;
SAC.hx2 = hx2;
SAC.hy2 = yh2;
SAC.mr1 = mr1;
SAC.mr2 = mr2;

% gaussian fit results:
SAC.Gfit1 = fp11; % parameters
SAC.niter1 = niter11; % iterations
SAC.chisq1 = chisq11; % fit error
SAC.err1 = err11; % parameter estimate error
SAC.dep1 = dep11; % dependency between parameters
SAC.fwhm1 = fwhm1;
SAC.NPH1 = max(yg1);
SAC.binw1 = binw;
%
SAC.Gfit2 = fp12;
SAC.niter2 = niter12;
SAC.chisq2 = chisq12;
SAC.err2 = err12;
SAC.dep2 = dep12;
SAC.fwhm2 = fwhm2;
SAC.NPH2 = max(yg2);
SAC.binw2 = binw2;

%CONTROL(sf).SAC = SAC; % save in the database.

figure(hf);
clf;


% text area
subplot('position', [0.1, 0.90, 0.8, 0.095]);
axis([0,1,0,1])
axis('off')
%ht(1)=text(0,0.80,sprintf('%-12s R[%d:%d]     %-8s  [%s]',DFILE.filename, DFILE.frec, DFILE.lrec, CONTROL(sf).protocol, date), 'Fontsize', 10);
%set(ht(1), 'interpreter', 'none'); % un-TeX the line - this is a filename and won't have tex chars, but might have an underscore.
%   text(0,0.6,sprintf('Solution:%-12s  gain:%4.1f  LPF:%4.1f kHz', CONTROL(sf).solution, DFILE.igain, DFILE.low_pass(1)), 'FontSize', 8);
%   text(0,0.4,sprintf('Ihold:%6.2f %s    RMP: %6.2f %s, Rin: %8.3f M\\Omega', ...
%      CONTROL(sf).iHold,CONTROL(sf).I_Unit, CONTROL(sf).Rmp, CONTROL(sf).V_Unit, CONTROL(sf).Rin), 'FontSize', 8);
   text(0,0.200,sprintf('Window1: %.1f-%.1f  mean rate: %.2f s/s s/s', start1, start1+dur1, mr1), 'FontSize', 8);
   text(0, 0.000, sprintf('G1: A=%.2f (%.2f) S=%.3f (%.3f)  NPH = %.2f  FWHM1 = %.3f', ...
       fp11(2), err11(2), fp11(4), err11(4), SAC.NPH1, SAC.fwhm1), 'Fontsize', 8);
   text(0, -0.200, sprintf('G2: A=%.2f (%.2f) S=%.3f (%.3f) NPH2 = %.2f  FWHM2 = %.3f', ...
       fp12(2), err12(2), fp12(4), err12(4), SAC.NPH2, SAC.fwhm2), 'Fontsize', 8);
   

subplot('position', [0.1, 0.075, 0.8, 0.320]);
%bar(sqrt(hx1), yh1, 'histc');
semilogx(hx1c, yh1, 'ks', 'MarkerFaceColor', 'k', 'MarkerSize', 3.5);
hold on
semilogx(hx1c, yg1, 'r');
%semilogx(hx1c, fitc11, 'g');

u1 = get(gca, 'Ylim');
h1 = gca;
set(gca, 'Xlim', [0 100]);
xlabel('Delay (ms)');
ylabel('Normalized Peak Height');


subplot('position', [0.1, 0.480, 0.8, 0.320]);

%bar(sqrt(hx2), yh2, 'histc');
semilogx(hx2c, yh2, 'ks', 'MarkerFaceColor', 'k', 'MarkerSize', 3.5);
hold on
semilogx(hx2c, yg2, 'r');
semilogx(hx2c, fitc12, 'g');

u2=get(gca, 'YLim');
h2 = gca;
set(gca, 'Xlim', [0 100]);
xlabel('Delay (ms)');
ylabel('Normalized Peak Height');
if(u1(2) > u2(2))
    u2(2) = u1(2);
end;
u2(1) = 0;
set([h1 h2], 'Ylim', u2);
y0 = u2(2);
x0 = 100;
axes(h1);
text(x0, y0, sprintf('Window 1 (%.1f - %.1f ms)', start1, start1+dur1), ...
    'horizontalalignment', 'right', 'verticalalignment', 'top', 'fontsize', 9);
box off
%axes(h2);
%text(x0, y0, sprintf('Window 2 (%.1f - %.1f ms)', start2, start2+dur2), ...
%    'horizontalalignment', 'right', 'verticalalignment', 'top', 'fontsize', 9);
%box off

if(pflag)
    orient landscape
    print;
end;
% write out text files to read in elsewhere to plot the data cleanly.
if(strcmp(CONTROL(sf(1)).E_C, 'Flat'))
    ht = fopen('/users/pmanis/desktop/popsacf1.txt', 'w');
    fprintf(ht, 'tf,df1,gf1\n')
else
    ht = fopen('/users/pmanis/desktop/popsacn1.txt', 'w');
    fprintf(ht, 'tn,dn1,gn1\n')
end;
for i = 1:length(hx1c)
    fprintf(ht, '%f,%f,%f\n', hx1c(i), yh1(i), yg1(i));
end;
fclose(ht);
if(strcmp(CONTROL(sf(1)).E_C, 'Flat'))
    ht = fopen('/users/pmanis/desktop/popsacfa.txt', 'w');
    fprintf(ht, 'tf2,df2,gf2\n')
else
    ht = fopen('/users/pmanis/desktop/popsacna.txt', 'w');
    fprintf(ht, 'tn2,dn2,gf2\n')
end;
for i = 1:length(hx2c)
    fprintf(ht, '%f,%f,%f\n', hx2c(i), yh2(i), yg2(i));
end;
fclose(ht);

function [y] = gaussfunc(x, fp)
%
% calculate a gaussian based on x and FP
%
y = fp(1) + (fp(2)/(fp(4)*sqrt(2*pi)))*exp(-((x-fp(3)).^2)/(2*fp(4)^2));


function printrasterplotfile(filename, index, Y)
hr = fopen(sprintf('/users/pmanis/desktop/%s', filename), 'w');
ml = 0;
for i = 1:length(Y)
    fprintf(hr,'tr%ds%d,rn%ds%d,',i, index,i, index);
    a = length(Y{i});
    if(a > ml)
        ml = a;
    end;
end;
fprintf(hr,'\n');
for j = 1:ml
    for i = 1:length(Y)
        if(j > length(Y{i}))
            fprintf(hr,',%d,',i);
        else
            fprintf(hr, '%f,%d,', Y{i}(j),i);
        end;
    end;
    fprintf(hr, '\n');
end;
fclose(hr);
